Skyrmion in spinor condensates and its stability in trap potentials 
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A necessary condition for the existence of a skyrmion in two-component Bose-Einstein condensates 
with SU(2) symmetry was recently provided by two of the authors [Phys. Rev. Lett. 97, 080403 
(2006)], by mapping the problem to a classical particle in a potential subject to time- dependent 
dissipation. Here we further elaborate this approach. For two classes of models, we demonstrate 
the existence of the critical dissipation strength above which the skyrmion solution does not exist. 
Furthermore, we discuss the local stability of the skyrmion solution by considering the second-order 
variation. A sufficient condition for the local stability is given in terms of the ground-state energy of 
a one- dimensional quantum-mechanical Hamiltonian. This condition requires a minimum number 
of bosons, for a certain class of the trap potential. In the optimal case, the minimum number of 
bosons can be as small ~ 10*. 
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I. INTRODUCTION 

Topological defects often play a fundamental role in 
our understanding of phases of matter and the tran- 
sitions between them. The best understood examples 
are probably vortices and vortex-loops in superfluids and 
0(2) magnets in two and three dimensions, which are re- 
sponsible for the very existence of the high-temperature 
phase, and completely determine the universality class 
of the phase transition [l|. The next in order of com- 
plexity are the topological defects in the 0(3)-symmetric 
Heisenberg model, which allows skyrmions in two dimen- 
sions |2] and hedgehogs-like configurations (a point de- 
fect around which spins point outward) in three dimen- 
sions 0]. While the former only renormalize the cou- 
pling constant, the role of the latter is less clear [J. All 
of the above, however, represent configurations topolog- 
ically distinct from vacuum, which provides them with 
local stability. 

In this paper, we study the stability of the topo- 
logically non-trivial skyrmion configuration in theories 
with 0(4) symmetry. Such a symmetry arises in com- 
plex condensates with an internal spin-l/2-like quantum 
number Q, for example. Realizations of such spinor 
condensates are found in models of inflatory cosmol- 
ogy [1], Bose-Einstein condensation of ^''Rb Q, and 
bosonic ferromagnetism [s!, 's'l , and in effective theories of 
high-temperature superconductivity [lo| and of decon- 
fined criticality [11]. The Higgs sector of the Weinberg- 
Salam model of electroweak interactions represents an- 
other closely related example, with a spinor condensate 
coupled to gauge fields. In our problem in three di- 
mensions, there exists a topologically non-trivial map- 
ping to the order-parameter space, thanks to the fact 
that the third homotopy group of is the group of inte- 
gers. However, the topology alone turns out to be insuf- 
ficient to guarantee local stability of the skyrmion. This 
may be understood already in terms of the classic Derick 



theorem [1^. Recently, a more general proof that the 
skyrmion cannot be a stationary point of the action for 
the spinor Bose-Einstein condensate (BEG) in free space 
was given by two of the present authors [13], based on 
the analogy between the Eulcr-Lagrange equations and 
the classical mechanics of a particle in a time-dependent 
dissipative environment .14] . The advantage of this alter- 
native point of view at the old problem is that it provides 
one with a simple way of constructing the external poten- 
tials which would indeed lead to skyrmion as the solution 
of the Euler-Lagrange equations. In Ref. [l^, three such 
special potentials were presented. In this paper, we fur- 
ther develop this approach. 

First we analyze the construction of the solution based 
on the ansatz proposed in Ref. [l^. There, a time- 
dependent dissipation which is odd in time was intro- 
duced, to allow an odd solution. However, we find that 
the symmetry argument does not always work, and there 
is a critical dissipation strength above which the odd so- 
lution no longer exists even for an odd dissipation. Next, 
we study the stability of the skyrmion in a generalized 
class of the potentials introduced in [l^ with respect to 
small variations. We map the problem to an effective 
quantum-mechanical eigenvalue problem and determine 
the region of local stability in the parameter space. A 
particularly interesting result of our analysis is that the 
stable skyrmion requires a minimal number of particles 
in the trap, estimated here to be ^ 10* in the optimal 
case. 

The paper is organized as follows. In SecHIl we present 
basic formulation of the problem. A classical equation of 
motion with time-dependent dissipation determines the 
skyrmion solution. On the other hand, a quantum me- 
chanical eigenvalue problem determines the local stability 
of the solution. In Sec. lIIIl we discuss the construction of 
the skyrmion solution based on the odd- function ansatz. 
We determine the critical dissipation strength, which sep- 
arates the region with and without a skyrmion solution. 
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In Sec. IIVI we present the numerical solution of the equa- 
tion of motion for a few cases. The existence of the crit- 
ical dissipation strength, as well as related theoretical 
predictions, is confirmed numerically. Furthermore, the 
local stability of the obtained solution is also analyzed 
by solving the quantum-mechanical eigenvalue problem. 
The minimal number of bosons required to satisfy the 
sufficient condition for the local stability is numerically 
obtained as a function of parameter. Section|V]is devoted 
to summary of the paper. 



II. GENERAL DISCUSSION ON THE 
SKYRMION SOLUTION 

A. Basic equations 



less parameters as 
S= I dT i(fr C 



£ = *t I 5^ _ 1 V2 - 1 + V{f)\ * + ^ {¥-^f . (5) 



Let us focus on classical field configurations indepen- 
dent of the imaginary time. It will prove convenient to 
represent the r- independent field ^ df) by an amplitude 
/(f) and a two-component complex spinor configuration 
a(f): 5'c('r) ~ f{f)a{f). The spinor a(f ) is normalized 
as adf)a{f) = 1. In terms of / and a, we can rewrite the 
effective action for the r-independent field configuration 



We begin by reviewing the derivation of the skyrmion 
solution for the two-component (spinor) BEC in an ex- 
ternal potential, formulated previously in Ref. The 
derivation is based on the Euler-Lagrange equations, and 
leads to the necessary condition for the existence of the 
skyrmion solution. We will then proceed to examine 
the stability of the skyrmion by taking into account the 
second-order variation. 

Let us consider the two-component bosons in three- 
dimensional continuum space in the external confinement 
potential V{r). The system can be described by the fol- 
lowing effective action via the path integral formalism [l| : 



2m 



V"^ - H + V{r) $ 



(1) 



where $^(t, r) = {^\{t , r) , ^'^{'^ : ''')) is a two-component 
bosonic field with the mass m, which satisfies the periodic 
boundary condition $(0, r) = <&(/?, r), where /3 = 1/k-QT. 
The bosons interact with each other by the repulsive con- 
tact interaction [/ > 0. is the chemical potential. In 
addition to the usual U(l) symmetry corresponding to 
the conservation of the number of the bosons, the system 
is invariant with respect to the global SU(2) transforma- 
tion of the boson field <& $'(t, r) — U^{T,r), where 
U is an SU(2) matrix. As usual, the trap potential V{r) 
has been included in the definition of the action or the 
corresponding Hamiltonian. 

Next, we introduce the dimensionless parameters and 
fields by rescaling: 



y(f ) = /i- V(r V), 



(2) 
(3) 
(4) 



where ^ = h/y/WiJI is the constant length scale. Effective 
action ([T]) may be written now in terms of the dimension- 



(6) 



where we have used the relation V(a'a) = (Va')a + 
a^(Va) = deduced from the normalization condition 
for the spinor a(f). The stationary state also needs to 
satisfy the boundary condition 



lim \f\'^f{rf\/a{f) = 0. 

|r I — >-oo 



(7) 



It guarantees the stability of the solution with respect to 
small rotations of the spinor a(f ) at the infinitely remote 
boundary of the system. 

Let us take the variation in action ^ around the clas- 
sical field ^'c(r) using the following expressions for the 
amplitude and the spinor configuration: 



fir) 
a{f) 



: /o(f ) + <5/(f ), 
ao(f) -I- (5a(f), 



(8a) 
(8b) 



where Sf{f) and Sa{f) are the variations around 5'c(f ) 
for the density profile and the spinor configuration, re- 
spectively. /o(^') and ao(T") are the amplitude and the 
spinor of some stationary field configuration. In par- 
ticular, ao(f) will assume a form corresponding to the 
skyrmion solution, which is to be defined shortly. Sub- 
stituting into the action, action ^ can be written as 



(9) 



where Ci is the Lagrangian density related to the ith or- 
der of the variation 5f{r) and 5a{r). The Lagrangian 
densities up to the second-order variation are then ex- 
pressed as 
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^0 = ^ ( V/o ) + ^/o'v4 • Vao + {V- + 



Ci = 

C2=Sf 



V^/o + {v4 • Vao + 21/ - 2} /o + 2C7/3] Sf 

5f + 



1/ - 1 + 3^7/2 - • Vao 



^ L (/o "^4) + /i-c 
/o (V,5at) + 2 (Vaj,) (5/] • [/o (vJa 



(10a) 
(10b) 



(Vao) Sf 



(10c) 



r 



B. Mapping to a problem in classical mechanics 

Setting J(Pf£i — for any Sf{f) and 6a{f) leads to 
the Euler-Lagrange equations giving the extremum of the 
action: 



ivaJj-Vao + V^-l 



/o + 2C//3 = (11) 



2V • [/o^Vao 
= "4V. (/o^Vao) +{V. (/o^Vaj)} 



ao 



(12) 



The latter equation follows by recalling that a^a = 1, so 
that aj(5a + Sa^ao = 0. 

For simplicity, we further assume that the external 
potential is spherically symmetric, V{r) = V{r), where 
f = |r|. Then, in terms of the spinor configuration, we 
can adopt the most general ansatz with the same spher- 
ical symmetry 



/o(r) = /o(r), 

/ sina;(f)cos0 — icosa;(f) 
I sinti>(f) sin 0e^"^ 



ao{f) 



(13) 
(14) 



where r — {f sin 6 cos (pjf sin 9 sin (j),f cos 9). At the in- 
finitely remote boundary, we impose u^r = 00) — Nn, 
where N is any integer. If we additionally adopt the 
boundary condition w(f = 0) = 0, and uj{f) changes 
from to Ntt as r changes from to infinity, ansatz (fT^ 
means that the spinor configurations wraps the three- 
dimensional sphere 5'^ N times, and thus corresponds to 
the skyrmion solution. Hereafter, we restrict the discus- 
sion to the simplest case of = 1. 

Let us then consider Euler-Lagrange equations pT|) 
and in]) for iV = 1 skyrmion and (dH). Substituting 
Eqs. and HH) into Euler-Lagrange equations (fTTj) 
and l|12p. we obtain the following differential equations: 



ffo 



+ 2 



df"^ 



24/0 
f df 

dhj\ ' 
'dfj 



^dfo 
fo df 



dio 
df 



sin 2uj 



V-l fo + 2Uf^^0. 

(15) 

= (16) 



To analyze these differential equations, it is very conve- 
nient to change the variable from f to t = In f. Then, as 
0<f<oo , —oo<t<oo. Rewriting Eq. (|16p in terms 
of t, we obtain 



dW , ^dui 



(17) 



a uj 

where u){t) = w(e*) = Lo{f). This equation can be re- 
garded as describing classical motion of a particle in an 
external potential with dissipation. W(ijj) and rj{t)^ re- 
spectively, correspond to the potential energy and the 
"time" -dependent dissipation, given by 

icos2w, (18) 



W{u)) 



1 



dt 



In/o', 



(19) 



where /o(i) = /o(e') = fo{f)- In terms of = 1 
skyrmion solution, boundary conditions for w(f) may be 
written as i^(0) — duj/df{Q) = at the origin of the 
three-dimensional space, and u;(oo) — n, duj/df — a,t 
the infinitely remote boundary. These boundary condi- 
tions in terms of t translate into 



diu 



(20) 



Integrating equation of motion (fT7|) with respect to t, 
and imposing the set of boundary conditions (j20p . one ob- 
tains the following necessary condition for the existence 
of the skyrmion solution: 



dt T]{t) 



duj\' 
It J 



0. 



(21) 



The condition implies that the total integrated dissipa- 
tion in the problem vanishes [lB|- It is now clear that 
a skyrmion solution exists only when the density profile 
takes a special form, so that the solution Lu{t) satisfies 
condition pT|) . 



C. Stability of the skyrmion and a 
quantum-mechanical eigenvalue problem 

The integral condition provides only the necessary con- 
dition for the existence of the skyrmion solution, but does 
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not guarantee its stability. Here, we analyze the second- 
order variation, to obtain a further condition for the sta- 
ble skyrmion solutions. 

Let us assume that an appropriate external potential 
V{f) is given so that the Euler-Lagrange equations allow 
a skyrmion solution fo{f) and ao(f). For the obtained 
skyrmion solution to be stable against local variations, 
the second-order variation in Lagrangian (jlOcp has to be 
positive. Obviously, the second term in Eq. (|10c[) is al- 
ways positive. The positivity of the first term in Eq. (jlOcp 
for arbitrary 5f{f) and 6a{f) is thus a sufficient condi- 
tion for the local stability of the skyrmion. The first term 
being a quadratic form of Sf, its positivity for an arbi- 
trary Sf is equivalent to positive definiteness of the linear 
operator 

H = -^W^ + VMf), (22) 

where 

V.sir) = V{f) - 1 + 3C//2(f) - ^Val ■ Vao- (23) 

Positive definiteness of H means that all the eigenvalues 
of H are positive. In fact, H may be interpreted as the 
quantum-mechanical Hamiltonian for a single particle in 
the external potential V^ff . The (sufficient) condition for 
the stability then corresponds to the ground-state energy 
of the Hamiltonian being positive. 

The external potential V{f) depends on both the 
spinor part and the density profile of the skyrmion, as set 
by the Euler-Lagrange equations. According to Eq. pT|) . 
Vc([{f) is also spherically symmetric. 

We should mention that even if the ground-state en- 
ergy of Hamiltonian is negative, it is possible that 
second-order variation (|10cP is still positive, if the posi- 
tive second term is large enough. We, however, will be 
unable to say more about this issue, and our discussion 
will be limited to the sufficient condition for the skyrmion 
stability formulated above. 

III. CONSTRUCTION OF SKYRMION 
SOLUTIONS FOR TRAPPED BECS 

A. Odd-function ansatz 

Let us discuss a few concrete examples of skyrmion 
solutions. In a usual formulation, we seek a solution 
for a given trap potential V{r). However, as discussed 
in Ref. ^13], a generic trap potential does not allow a 
skyrmion solution. Thus, we solve the problem backward: 
first we determine the dissipation ri{t) so that equation 
of motion p7|) has a solution uj{t) which represents a 
skyrmion. The assumed dissipation determines the den- 
sity profile of bosons. Finally, the trap potential V{r) is 
determined so that it reproduces the chosen ri{t). 

For convenience, here we introduce the new variable 

y^cs-f- (24) 



Equation of motion (|17|) is then rewritten as 

g = -.(t)f -sin2,. (25) 
Boundary conditions (pn|) read, in terms of y, 

y(-oo)--|, ^(-^) = 0, (26) 

y(oo) = +|, ^M = 0- (27) 

We observe that the boundary conditions at i = -l-oo 
would be automatically met if y obeys boundary condi- 
tions (j26p at t = — oo and y{t) is an odd function: 

y{t) = -y{-t). (28) 

This is sufficient to satisfy the original conditions ([^5]) 
and (|27p . but not necessary. However, here we focus on 
finding the solutions y(t) which are odd in t, because it 
is easier than solving the general problem. 
Equation (|25p is invariant under 

t ^ -t, y ^ -y, (29) 

if ri{t) is also odd in t. Thus we choose an odd ri{t) so 
that we can find an odd solution y{t). Although this 
is a somewhat restrictive choice, it is a useful ansatz in 
construction of skyrmion solutions. 

However, it turns out that some choice of odd r]{t) ac- 
tually does not allow an odd solution y{t) which satisfies 
the boundary conditions at i = — oo. Roughly speaking, 
if the dissipation is too strong, the particle which starts 
at y — — 7r/2 with a vanishing speed at t = — oo can not 
reach y = at t — 0. We will examine two forms of odd 
rjit) as examples: 

T^{t) =n[e{-t - 1) - 9{t - 1)], (30) 
ry(<) = — rttanht, (31) 

where n is a positive parameter. As we will demonstrate, 
for each case there is a critical parameter ric; an appro- 
priate odd solution exists only if n < ric- 

The existence of an odd solution can be discussed in 
the following manner. First we impose boundary con- 
ditions (pS]) only At t = —oo. Of course, they do not 
completely fix a solution but allow a family of different 
solutions. This is evident by recalling that the trivial 
solution y{t) — — 7r/2 satisfies Eq. ([26)) . 

Then we attempt to find, among the solutions, the one 
which satisfies y{0) = 0. This would be the desired odd 
solution. The particle falls off the hill and approaches 
the potential minimum y — 0. If the particle is still 
rolling down the hill (— 7r/2 < y < 0) when t — 0, the 
particle is only accelerated by the negative friction for 
t > and it must eventually reach y = at some pos- 
itive t. Namely, for any dissipation strength, there is a 
solution which satisfies y{t) = at a positive t. On the 
other hand, while y < 0, the velocity dy/dt should al- 
ways be positive. Thus the smallest solution i of y(t) = 



5 



changes continuously. Therefore, if there is another solu- 
tion which satisfies y{t) = at a negative t, there must be 
an odd solution with y(0) =0 thanks to the intermediate 
value theorem. 



B. Constant- r; regime 

In either case of Eq. ([50]) or (PT|) . we observe that, for 
t < -1, 



7?(i) - n. 



(32) 



This simply represents the constant-dissipation coeffi- 
cient. The equation of motion 



d y 



dy 

-n sm 2y. 

dt 



(33) 



in this regime is independent of time. Thus, for any solu- 
tion y(t), the translated solution y{t+T) is also a solution 
for any r. Let us discuss the solutions of this equation. 
This will turn out to be useful in determining the critical 
parameter ric for the original equation of motion (j25p . 

As we will see later, when n ric, the particle comes 
very close to the minimum (|y| <^ 1) while still in the 
constant-ry regime {t <ti —1)- For small y, we may use 
the linearized equation of motion 



d'^y dy 



(34) 



instead of the full nonlinear equation p3|) . Its solution 
can be easily obtained as 



y{t) = Aie-^i(*+^) -I- A2e-^'^'+^\ 



where 



Ai = 



A, 



n — \J V? 



(35) 

(36) 
(37) 



One can choose an arbitrary large positive t, thanks to 
the translation invariance in time. On the other hand, 
taking a large negative r (for a fixed t) makes y large, 
and may invalidate linear approximation ()34|) . 

When n < 2^/2, Ai,2 are a complex conjugate pair and 
the solution represents a damped harmonic oscillation. 
In this case, there is obviously a solution which reaches 
y = within Eq. ((33)) . Namely, there is a solution which 
satisfies y{t) = at a negative t. As we discussed in 
Sec. nil AI the intermediate value theorem assures that 
there is an odd solution with y{0) = in this case. This 
means that ric > 2\f2. 

Thus, in the following, we focus on n > 2\/2. Then, the 



first term 



is the leading one in the r 



The second term e"'*'^'*^'^'' vanishes more quickly, but the 
subleading contribution determines the critical exponent 
as we will show later. 



In fact, in discussing the subleading contribution, we 
must also consider the nonlinear effects which were ig- 
nored in Eq. ([34| . We introduce the scaling of y by the 
replacement y — > ay. We are interested in the limit in 
which the original y is small, namely a — » 0. The equa- 
tion of motion now reads 



d_y_ 

dt^ 



dy 
I— 
dt 



1 



sin (2ay). 



(38) 



Considering the limit a ~* and retaining only the lead- 
ing nonlinear term, we obtain 



ay 

dt^ 



dy 
I— 
dt 



2?/ 



(39) 



We consider a series expansion of the solution y in terms 
of a, which can be regarded as a perturbative expansion 
of nonlinear effects. 

The lowest order j/^'^' is given by the solution of lin- 
earized equation ((34|) . The next order y*^^) is of O(a^), 
and is given by a solution of 



^2 (1) ^ (1) 4 



(40) 



This is an inhomogeneous linear differential equation on 
y^^^ for a given y^^\ which can be solved by a standard 
method. Taking the as general solution ([55)1 of the 
linear equation, we find the special solution 



2Ai' 



3Ai(A2 — 3Ai) 



-3Ai (t+T) 



(41) 



where only the leading term is given. The general solu- 
tion also contains solutions of the corresponding homo- 
geneous equation. However, they have the same form as 
y^^^ and can be ignored in the following. 

Combining with the solution of the linearized version 
[Eq. (|35p] , the leading and next-leading terms in the r 
oo limit can be written as 



where 



yit) = Aie-^i(*+")(l + Ce-'^^(*+")), (42) 



AA EE min(A2 - Ai,2Ai). (43) 



Namely, when A2 > 3Ai, the next-leading term comes 
from the nonlinear effect instead of the term proportional 
to e--^2(t+r)^ 

Recalling that the particle approaches y = from 
y = — 7r/2, Ai < 0. The constant C is determined by 
the solution of nonlinear equation of motion p3p with 
the constant dissipation and boundary conditions (|26p . 
The solution gives the effective initial conditions for the 
linearized equation. 

Numerically solving Eq. ([55)) with Eq. (|^^ . we find 
that, the ratio —{dy/dt)/y increases monotonically, as 
shown in Fig.)!] when n > 2^2. This observation implies 
that C < in linear regime (|32]). 
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FIG. 1: (Color online) The ratio —(dy/dt)/y as a function 
of time t, within constant-r; equation (I33|l and with boundary 
conditions (|26p . In the actual numerical calculation, a small 
initial velocity is given instead of zero in Eq. (|26p to avoid 
getting only the trivial solution. The result is robust against 
the change in the small initial velocity. In the figure, we 
show the result for = 4 as example. The ratio —{dy/dt)/y 
increases monotonically, asymptotically approaching Ai = 2 — 
V2 ~ 0.585786. We obtained similar results for other 77 > 



In fact, when the nonlinear effect gives the next-leading 
term, it follows from Eq. (PT|) that 



C = ^c 



< 0. 



(44) 



3Ai(A2 -3Ai) 

When the subleading contribution within the linear equa- 
tion gives the next-leading term, we do not have a proof 
but C < seems certain from numerical results. 

Furthermore, the numerical results indicate that C is 
sufficiently small so that y never reaches 0, when asymp- 
totic expression is valid. This, of course, does not 
mean that there is no solution satisfying y(0) = in the 
original Eq. (j25p . in which the dissipation is turned off 
around t ^ 0. 



C. Critical parameters in the step- function case 

Let us consider the step- function case of Eq. ([50)1 . 
Here, the solution of Eq. discussed in Sec. IIII Bl gives 
the "initial condition" at t — ~1 for the equation with- 
out the dissipation. If the particle is very close to the 
minimum {\y\ <^ 1) at t =^ —1, the successive motion is 
just a harmonic oscillation with the angular frequency of 
^/2. We will show later that for n '--^ ric, ly] <C 1 indeed 
holds at t = —1. 

The phase C of the oscillation is given as 



tanC = y\/2 ( ^ 



(45) 



When n — > ric, the particle just manages to reach y — 
at t = 0. For that, we need to give the optimal initial con- 
dition at t = — 1, namely, —(dy/dt)/y with the maximum 



possible value. As discussed in Sec. IIII Bl ~-{dy/dt)/y for 
Eq. (j33p monotonically increases. Thus the optimal ini- 
tial condition [maximum possible —{dy/dt)/y\ is realized 
by letting the particle spend infinite time around y ^ Q 
before t = —1, namely, by taking t ^ 00. In this limit, 
the first term in Eq. (H^ dominates and the initial con- 
dition at t = — 1 is given by 



tan Co = — 



V2 



(46) 



We emphasize that T depends only on the ratio between 
the "initial" velocity and "initial" coordinate on i = — 1, 
which converges to a finite value in the limit t ^ 00. 
The time required to reach the minimum {y ~ 0) in the 
harmonic oscillation is 



T^Ml tan- ^. 



(47) 



The critical dissipation coefficient ric in this problem 
is thus given by 



T = ^tan^i ^ 1, 
\/2 Ai 



(48) 



so that the particle arrives at y = at t = 0. Therefore 
we find 



2V2 
in2V2 



9.18107. 



(49) 



In the limit n — )■ iic, t — > 00, and thus both ?/(— 1) 
and dy/dt{~l) vanish. This implies that the velocity 
of the particle when it reaches the potential minimum, 
dy/dt{0), also vanishes. Let us discuss its critical be- 
havior, namely, how dy/dt(Q) depends on ric — n when 
n < fic- 

For n < 71c, if we take t — > cxd the particle reaches 
the minimum before t = because T < 1. By letting the 
particle spend less time around y ~ by taking a smaller 
T, we can change the initial condition at t = — 1 so that 
Co is smaller than the optimal value Therefore, for 

n < ric, there is a solution which reaches the minimum 
at i = 0. Combining 



f ^ I (t = -1) = tan^ 



(50) 



with Eq. (im, for small Uc ~ n wc find 

(n-n,) + (A2-Ai)Ce^'^^(^~i)-0. (51) 
This implies that 



dXi 
dn 



e '^'^'^ oc Uc — n. 



(52) 



As a consequence, 
dy , 



y{-l) cx ^(-1) cx e-^^- (x (n, - n)^^/^\ (53) 
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which leads to 



^(0)(xK-n)^^/^\ 
at 



(54) 



In the present case, for n ~ ric, A2 > 3Ai, namely, AA 
2Ai. Thus we obtain 



dy 
dt 



^(0)(xK-n)i/^ 



(55) 



For n > Tie, even under the optimal condition r ^ cx3, 
the particle cannot reach y = at t = 0. In this case, 
the odd solution does not exist. 



D. Critical parameters in the tanh case 



Now let us consider tanh case (|3T|) . The mathematics 
is somewhat more complicated but the physics is quite 
similar to the previous one. 

As in the previous problem, for n ~ rtc, we can assume 
that the particle comes very close to the minimum y ^ Q 
at negative time. The linearized equation of motion, with 
the full time dependence of the dissipation, reads 



dy 

-n — tanht — 2y. 
dt 



(56) 



This equation has the general solution 

y =CiP^;ff^^^\tanht)il - tanh^ i)-"/^ 

+ C^Q^n^f^^^^ (tanh t){l- tanh^ t)-"/4, (57) 

where Ci^2 are constants, and P^i"^^ and Qji'^ are associ- 
ated Legendre functions. 

Let us first discuss the asymptotic behavior of the 
above solution in the limit t —00. In this limit, equa- 
tion of motion reduces to Eq. and thus the 
solution should be equivalent to Eq. (|35p . In fact, using 



tanht l + 2p/* + 0(e'^'). 



(58) 



and the asymptotic expansion of the associated Legendre 
functions, we have confirmed that Eq. (|57p coincides with 
Eq. (|35p in the limit t —^ —00. 

Now, the similar discussion as in Sec. lIII Cl applies here. 
Namely, the numerical solution implies that —{dy/dt)/y 
increases monotonically. Thus, the optimal condition for 
reaching ?/ = is realized when the particle spends infi- 
nite time around y ~ before t ^ Q. This can be done by 
replacing t with t + t and taking r — > 00. In this limit, 
the asymptotic behavior of the solution in the t <C — 1 
regime is dominated by e~'^i'^*+'^^ 

In the following, we demonstrate that 



3. 



(59) 



To show that, let us set n 
solutions reduce to 



3. The two independent 



P3'J/^)(tanht)(l - tanh^ t)-3/4 
g^/f ^(tanht)(l - tanh^ t)-^/^ 



1 



2V27r 



Alt 



(60) 



(61) 



-00. This implies that, under the op- 

(1/2) 



in the limit t 

timal condition, the solution consists only of the 
term. 

On the other hand, for n = 3, the Taylor expansion 
around i = reads 

Pg^Jf ^(tanhi)(l - tanh2<)-3/4 _ _ ^^o{t^), (62) 

Q^y/^ (tanh t)(l - tanh^ t)-'^''^ ~ - yphit + 0{t^). (63) 

Namely, the solution consisting only of the Q term just 
crosses y — Q on t — Q. By perturbing the solutions 
around n = 3, it can be shown that the solution with 
y(0) — exists for n < 3. This means that 71 = 3 is 
indeed the critical value Uc. 

As in the case in Sec. IIII CI the velocity at the poten- 
tial minimum dy/dt{0) vanishes as n approaches ric from 
below. The critical behavior can be obtained by a sim- 
ilar argument, and is given by Eq. ([5i)) . In the present 
case, Ai = 1 and A2 = 2 < 3Ai. Thus we find the linear 
behavior 



^(0) cx (ric-n), 



(64) 



when n < 



E. Requirement of a finite number of bosons 

We consider the situation where all the bosons are con- 
fined in the finite space by external trap potential V{r). 
If it is to be realized in experiments, the total number 
of the bosons should be finite. This gives an additional 
requirement independent of the stability. 

In terms of the density profile /o^, the condition of the 
finite number of bosons is easily given as 



df /o^(r) = 47r / dr r^/o (r) < 00. 



(65) 



For the choice of the tanh dissipation [Eq. (|?T|) ]. the 
density profile can be obtained via Eq. (fT9)l as 



foir) = B 



(66) 



where _B is a positive constant which appears from the 
integration of Eq. (fT9|l . The density profile in the n = 1 
case was discussed in Ref. [ll]. r^/o (r) in the integral of 
Eq. ([65]) asymptotically behaves as ^ f-"+i for large f. 
Accordingly, for n < 2 the total number of the trapped 
bosons diverges, violating condition (j65p . Hereafter, we 
focus on the finite-bosons case, n > 2. 
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FIG. 2: (Color online) Trajectories of cij(t) in the case of 
r\(€) = — ntanht. The trajectories in f < are antisymmetric 
to that in t > 0, so we show only for t > 0. Here, 

calculating the trajectories numerically, Loit = 0) = 7r/2 is 
imposed as an initial condition. Now, although we define 
ri > 2, the trajectory for n = 2 is also shown. 




10"^ 10"^ 10"'' 

Uc — n 

FIG. 3: The velocity of the particle at t = as a function of 
ric — n, for the tanh dissipation [Eq. pip ]. The squares are 
numerical results obtained by the shooting method. The line 
is the best fit assuming a power law. The critical parameter 
Jic is determined by the analytic prediction. The good fit 
to the power law means that the analytic prediction of ric 
is consistent with the numerical calculation. Moreover, the 
exponent 0.979 ■ • ■ obtained by the fit is also consistent with 
the analytic prediction of unity. 



IV. NUMERICAL SOLUTION FOR A 
SKYRMION 

A. Solution of the equation of motion 

Equation of motion (jl7p determines (D(t) for a given 
choice of r\(t). While we have determined critical param- 
eters in Sec lIIIl unfortunately, the full solution cannot 
be obtained analytically. Thus, here we solve Eq. pT]) 
numerically. The numerical solution can be also used to 
check the analytical predictions on the critical parame- 
ters discussed in Sec. IIIII 

To reiterate, we seek a solution which satisfies the 
boundary conditions at t — ±oo [Eq. (|^D|) ]. For an odd 
T](t), which is the case we discuss in this paper, such a 
solution satisfies Eq. ((28)) . To find the odd solution for 
an odd ri{t), it is enough to require 

uiO) = 0, (67) 

together with either of the boundary conditions at < = 
— oo or i = oo in Eq. 

Based on this observation, we adopt the so-called 
"shooting method" in the numerical scheme, which is 
explained in Appendix [Xj In Fig. [2l we show the nu- 
merical result for the case of Eq. (PTjl . Here, we observe 
that the velocity at t = vanishes as the dissipation 
parameter n approaches ric — 3. This is indeed consis- 
tent with the analytic prediction on the critical param- 
eter Uc = 3 and on the critical behavior. To see this 
more clearly, in Fig. [31 we show the numerical result on 
du!/dt{0) = dy/dt{0) as a function oiric — n. The result is 
in good agreement with the analytic predictions. We also 
have made a similar comparison for step-function dissi- 
pation ([50)1 in Fig. [H and found agreement with analytic 
predictions ([^0)1 and ([55)1 as well. 




FIG. 4: The velocity of the particle at the potential minimum 
{t = 0) as a function of ric—n, for the step function dissipation. 
The line is the best fit assuming a power law to the numerical 
results, which are shown as squares. The critical parameter 
Tic is determined by the analytic prediction. The numerical 
result is again in good agreement with the predictions. 



B. Stability of the skyrmion and minimal number 
of bosons 

We have shown that the skyrmion as a solution of the 
Euler-Lagrange equation exists for n < Uc, and obtained 
the solution numerically. The spinor configuration uj(t) 
is given by the solution of a classical mechanics problem. 
The density profile fo{r) can then be obtained. 

The next step is to examine the stability of the ob- 
tained skyrmion. Note that there is a free parameter B 
appearing in density profile ([66)) . Substituting the ob- 
tained u){t) and fo{t) into Euler-Lagrange equation ([TB)l . 
the trap potential V{r) realizing the skyrmion is ob- 
tained. Then, we find that V(f) depends on U and B 



9 



400 



200 



4 



-200 



in 



in 



B = 
B = 100 
B = 500 
B = critical 



2000 



in 



2000 



-4000 



n = 2.0 
n = 2.5 
n = 2.7 
n = 2.9 



0.1 



0.2 



FIG. 5: (Color online) Plots of the eflFective potential VcS 
defined in Eq. (I23|) for given ri{t) = —ntanht. In the upper 
graph, V'off(?) with fixed n = 2.5 are plotted. In the lower 
graph, the n dependence of 14ft (?) with fixed S = is shown. 
For each value of n, there is a critical value of B which makes 
the ground-state energy of Hamiltonian (|22p exactly zero. In 
the region of B shown in the upper graph, the ground-state 
eigenvalues are negative. 




FIG. 6: The critical value Be as a function of n. The values of 
Be axe plotted on logarithmic scale. We find that Be increases 
with n. Eventually, in the limit n ^ 3, Be diverges. 






^10^ 



FIG. 7: The n dependence of the number of the trapped 
bosons m B = U~^Be for a given n. It has the minimal value 
of ~ 10"* around n = 2.25. 



as a function oi UB = B. This information on uj{r), 
/o(r), and V{f) determines the effective potential p3|. 
shown in Fig. O In order to obtain the stable skyrmion, 
B larger than the critical value Be is needed. Indeed, cal- 
culating the ground-state energy of quantum-mechanical 
Hamiltonian (p2|) by use of the numerical diagonalization 
method, it is found that as n increases, the critical value 
Be diverges as n approaches nc = 3, as shown in Fig [61 

Finally, let us investigate the total number of the 
trapped bosons for some values of B and n. Due to 
Eq. (|66p for the density profile, the number of bosons 
monotonically increases as a function of B: it is propor- 
tional to B. For a given n, the number of the trapped 
bosons is minimal at B = Be- The n dependence of the 
number of the bosons in ,B = ,Bc is then shown in Fig. [T] 
We find the minimum at n 2.25, with the value of 
~ 10*. It is comparable to the typical numbers in the 
experiments. 



V. SUMMARY 

In this paper, we have discussed the skyrmion config- 
urations of the two-component spinor BECs confined by 
a trap potential. The necessary condition for the exis- 
tence of the skyrmion has been formulated earlier by two 
of the authors [l3| . The Euler-Lagrange equation, which 
must be satisfied by the skyrmion, turned out to give 
an equation of motion for a fictitious classical particle 
subject to time-dependent dissipation. A skyrmion so- 
lution satisfying appropriate boundary conditions exists 
only for specially chosen trap potentials. Some of those 
trap potentials were constructed by considering dissipa- 
tion which is an odd function of time. It was expected to 
allow a solution of the equation of motion, which is odd 
in time. 

In this paper, we have developed this approach fur- 
ther. We have mainly considered two classes of solu- 
tions, obtained by taking the dissipation as proportional 
to the step function of time [Eq. ([50)1 ]. and to hyperbolic 
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f 

FIG. 8: The example of the trap potential V{f) reproduc- 
ing the skyrmion. The density profile under this potential is 
given by Eq.((66l) with n = 2.25 and B = 213.495. As seen in 
Fig. [71 for these n and B, the number of the trapped bosons 
is minimal (-^ 10^). 



tangent of time [Eq. ([3T|l ]. We have found that, even 
though these choices of dissipation are odd in time, they 
do not allow a skyrmion solution if the dissipation is too 
strong. We have obtained the critical parameter exactly, 
and also determined a critical behavior in the velocity at 
the potential minimum. These predictions are verified by 
numerical solution. 

Furthermore, we discussed the stability of the 
skyrmion, taking into account the second-order varia- 
tional theory. In our formulation the problem of the sta- 
bility of the skyrmion is mapped onto that of the sign 
of the lowest eigenvalue of a certain quantum-mechanical 
Hamiltonian determined by the skyrmion solution. 

For the tanh model, the density profile and the trap 
potential reproducing the skyrmion have several param- 
eters U , n and B, and we determined the region in the 
parameter space that leads to a stable skyrmion. For a 
given n, B — UB needs to be larger than a certain value 
for the skyrmion to be stable (see FiglB]). 

We have found that there exists a minimal number of 
bosons of ~ 10^ allowing a stable skyrmion solution. For 
illustration, the trap potential required to reproduce the 
skyrmion with this minimal number of particles is shown 
in Fig. [51 In order to connect our results to real experi- 
ments, it would be important to see the asymptotic form 
of the trap potential stabilizing the skyrmion. It can be 
analytically derived from asymptotic analysis of differ- 
ential equations (fTSl) and (fTHl) . In the generic case that 
the odd ri{t) asymptotically reaching =pn, respectively, as 
t — !■ ±cx), as seen in Appendix [Bl the form of the trap 
potential can be obtained as Eqs. (|B5[) and (|B6[) . respec- 
tively, for f <C 1 and f ^ 1. One can easily find that 
these equations are consistent with the result of the an- 
alytically solvable cases discussed in Ref. r]{t) = 
for n = and ri{t) = — tanht for n — 1. For n > 2, the 
leading terms in Eqs. (|B5p and (|B6p are both propor- 
tional to ^ 1/f^. Interestingly, the power of the leading 



terms is independent of n. In particular, in the n = 2.25 
case shown in Fig. [51 the analytically obtained asymp- 
totic form for f ^ 1, 1/f^, is in good agreement with the 
numerical result. On the other hand, in the case of the 
large-f side, because of the very large prefactor B ^ 10^, 
the subleading term remains effective. Thus, taking into 
account the leading term as well as the subleading term 
~ — S/f"+^, the analytical result for n = 2.25 matches 
the numerical one very well. 

There are many open problems which would deserve 
further study. In this paper we only discussed skyrmions 
in equilibrium. Dynamical aspects, such as relaxation to 
the skyrmion solution, are also important. The dynami- 
cal aspects would be even more crucial in possible phys- 
ical realizations of our proposal. As we have shown, the 
potential has to be fine-tuned to allow a stable skyrmion 
solution. In reality, it is impossible to construct a po- 
tential with infinite precision and thus there would be 
some error in the potential. This would make skyrmions 
absent, as a stable solution in equilibrium. However, we 
expect that, if the actual potential is close enough to the 
exact one with a skyrmion solution, the skyrmion would 
exist as a quasi-stable state with a certain lifetime. In or- 
der to discuss the feasibility of observation of a skyrmion, 
we would need to estimate the lifetime of the skyrmion 
in the actual potential. We hope progress will be made 
on these problems, and also on other directions related 
to our study. 
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APPENDIX A: SHOOTING METHOD 

In general, iterating the discretized time step based on 
the boundary condition, the equation of motion can be 
numerically solved. Based on this observation, we adopt 
the so-called shooting method in the numerical scheme 
as follows: 

1. Choose the initial velocity doj/dt{0) arbitrarily. 

2. Given the initial velocity and the initial position 
<D(0) ~ 0, solve equation of motion (fTTl toward 
t — oo. 

3. If the particle goes beyond the peak of the potential 
at a) = TT, the initial velocity was too large. 
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4. If the particle comes back without reaching the 
peak of the potential at Co — tt, the initial veloc- 
ity was too small. 

In a true solution, the particle should approach asymp- 
totically the peak of the potential uj = n as t oo. How- 
ever, the asymptotic behavior in t — > oo is quite sensitive 
to the initial velocity, and the numerical solution departs 
from the peak of the potential in either way, depending 
on the tiny difference of order of machine precision, in the 
initial velocity. The sensitivity is due to the instability 
of the particle at the potential peak, further enhanced by 
the negative dissipation coefficient. On the other hand, 
because of this sensitivity, we can easily obtain the cor- 
rect initial velocity in a high precision. 

In practice, the shooting method can be implemented 
as an efficient iteration using the bisection method as 
follows: 

1. Find the two values of the initial velocity {wf, }, 
so that Vi is "too small" and vY is "too large." 

2. On nth iteration, the correct initial velocity should 
be within the range {v^,v^). Thus, numerically 
solve the equation of motion with the midpoint ini- 
tial velocity = {v^ + v^)/2. 



3. If the midpoint is too small as the initial veloc- 
ity, set = v^', vY,+^ = vK, If v^" is too large, 
set instead v^^^ — v^, v^^i = v^/ . Go to step 2 as 
the (n -I- l)th iteration. 

In this way, the error in the initial velocity decreases 
proportionally to 2^" in the nth iteration. 



for t ^ — 1 , and 



-2tt 



duj 



(B2) 



for i 3> 1. The solutions satisfying the boundary condi- 
tions, lim(^_oo = and lmit^^^uj{t) = tt, can be 
easily obtained, respectively, as 



Die 



At 



D2e 



-At 



(i < -1 
it » 1) 



(B3) 



where A ~ 



> 0. Di and D2 are unknown 



constants. On the other hand, from the boundary con- 
ditions for r](t) and the definition of the dissipation 
ri{t) = 1 + ^ Iu/q, the asymptotic forms of the density 
profile for f ^ — 1 and t ^ 1 are easily obtained as 



Pit) oc 



g-(n+l)f ^ 



(B4) 



Taking asymptotic forms (|B3P and (jB4|l from the t 
representation to f representation hy t ^ Inf, and sub- 
stituting them into one of equations of motion (|15[) . the 
trap potential V{f ) can be derived. As a result, V(f ) for 
f <C 1 and f ^ 1 are investigated as 



V{P) = l-Dl- 



A^ + 2 

2f2-2A 



n2 - 1 U , , 

— ^-Si^, B5 



APPENDIX B: ASYMPTOTIC FORM OF TRAP 
POTENTIALS 

In order to drive the asymptotic form of trap poten- 
tials, we rewrite equation of motion p?|) for t <C — 1 and 
t ^ I. Then, in this odd-dissipation case r]{t) can be 
approximated as r]{t) ~ n for t —1 and r]{t) « —n 
for t ^ 1. In addition, as shown by the numerical re- 
sults, the solution would stand near Lu{t ^ —1) w and 
uj(t ^ 1) ~ TT. From these assumptions, equation of 
motion (|17p is linearized as 



d^w _ dto 



(Bl) 



for f <C 1, and 



V{r) = 1- D\ 



A^ + 2 

2^2+2A 



8f2 



U 



2 — 



l+n ■ 



(B6) 



for f 3> 1. Bl and B2 are unknown constants and are the 
prefactors of the density profile, respectively, for i ^ — 1 
and t 3> 1. In particular, in the case of ri{t) = — ntanht 
discussed in this paper, they turn out to be B appearing 
in Eq. ([55]) : Bi = B2 = B. The leading term is decided 
by the value of n. However, in Fig. [8l because of the 
very large prefactor BU ^ 10^, the subleading term still 
remains effective in shown region of f. 
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